#pragma once
#ifndef _RAYTRACER_OPTICS_H
#define _RAYTRACER_OPTICS_H

#include "numeric.h"
#include <maya/MFloatVector.h>


//---------------------------------------------------------------------------
// calculates the reflection vector of incident vector w/rt normal.
// incident vector points away from reflection point (so dot(I,N) > 0)
// incident and normal are assumed normalized
inline MFloatVector reflection(const MFloatVector& incident, const MFloatVector& normal)
{
  float norm_vec_dot = incident * normal;
  norm_vec_dot  = clampbot(norm_vec_dot, 0.0f);
  return 2*norm_vec_dot*normal - incident;
}

//---------------------------------------------------------------------------
// calculates transmission of incident vector w/rt normal, 
// from material with coeff. eta_i into one with coeff. eta_t.
// incident vector points away from reflection point (so dot(I,N) > 0)
// incident and normal are assumed normalized
inline MFloatVector transmission(const MFloatVector& incident, const MFloatVector& normal,
                                 float eta_i, float eta_t, bool& valid_result, float* ptheta_i=NULL, float* ptheta_t=NULL)
{
  // calculate theta_i
  float norm_vec_dot = incident * normal;
  float theta_i = acos(norm_vec_dot);
  if(ptheta_i){
    *ptheta_i = theta_i;
  }

  // calculate theta_t
  float ratio_eta_it = eta_i/eta_t;
  float sin_theta_t = sin(theta_i) * ratio_eta_it;
  float theta_t;
  if(sin_theta_t <= 1.0f){
    //sin_theta_t = clamptop(sin_theta_t, 0.99985f);
    theta_t = asin(sin_theta_t);
  } else {
    valid_result = false;
    if(ptheta_t){
      *ptheta_t = 0.0f;
    }
    return MFloatVector(0.0f, 0.0f, 0.0f);
  }
  valid_result = true;

  if(ptheta_t){
    *ptheta_t = theta_t;
  }

  // calculate T
  MFloatVector T1 = ratio_eta_it*(cos(theta_i)*normal - incident);
  MFloatVector T2 = cos(theta_t)*normal;
  MFloatVector T  = T1 - T2;
  return T;
  // there's another expression with just I.N and ratios and 1 sqrt

}

// should come out to same as transmission function above
inline MFloatVector transmission2(const MFloatVector& incident, const MFloatVector& normal,
                                 float eta_i, float eta_t)
{
  float ratio_eta_it = eta_i/eta_t;
  float i_dot_n = incident * normal;
  float inner_term1 = ratio_eta_it*i_dot_n;
  float inner_term2 = sqrt(1-((ratio_eta_it*ratio_eta_it)*(1-(i_dot_n*i_dot_n))));
  MFloatVector T = -ratio_eta_it*incident + (inner_term1 - inner_term2)*normal;
  return T;
}

inline float fresnel_kr0(float eta_i, float eta_t)
{
  return ((eta_i-eta_t)*(eta_i-eta_t))/((eta_i+eta_t)*(eta_i+eta_t));
}

inline float fresnel_kr(float kr0, float theta_i)
{
  float one_minus_cos_theta_i = (1-cos(theta_i));
  return kr0 + 
    ((1-kr0) * one_minus_cos_theta_i*one_minus_cos_theta_i*one_minus_cos_theta_i*one_minus_cos_theta_i*one_minus_cos_theta_i);
}

#endif
